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Abstract 

Event-driven molecular dynamics is a valuable tool in condensed and soft 
matter physics when particles can be modeled as hard objects or more gen- 
erally if their interaction potential can be modeled in a stepwise fashion. 
Hard spheres model has been indeed widely used both for computational 
and theoretical description of physical systems. Recently further develop- 
ments of computational techniques allow simulations of hard rigid objects of 
generic shape. In present paper we will present some optimizations for event- 
driven simulations that offered significant speedup over previous methods. In 
particular we will describe a generalization of well known linked list method 
and an improvement on nearest neighbor lists method recently proposed by 
us. 
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1. Introduction 

Systems composed of many particles can be modeled as hard rigid bodies 
(HRB) if excluded volume interactions are dominant and despite absence of 
any attraction they exhibit a rich phase diagram especially if their shape 
is non-spherical [Tj [2]. The spherical version of these models has already 
proven to be quite flexible and have been used to tackle, for example, several 
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biological problems [3H5]. The generalization to non isotropic object will 
increase even more their flexibility and applicability. 

Furthermore also attractive interactions between HRBs, provided that 
they are short-ranged and/or localized, can be properly modeled employing 
sticky spots (SS) [SHE]. 

Several numerical techniques have been developed in the past to perform 
molecular dynamics simulations of particles interacting with only excluded 
volume interactions. Dealing with hard bodies, the system is propagated in 
configuration space from one event to next, giving rise to so-called event- 
driven molecular dynamics (EDMD). The essence of these EDMD numerical 
algorithms involves the evaluation of the overlap between different objects 
[9HT2] or, equivalently, their geometrical distance [13] . 

In this paper we propose two optimizations for EDMD, namely: mul- 
tiple linked lists (MLL) method and nearest neighbor lists with null spots 
(SNL). MLLs method proved to be very useful in simulating mixtures of 
hard spheres with very different sizes, while SNL method offers a significant 
speedup for simulating particles with more complicated shapes, like hard- 
ellipsoids (HE)[lJ or superquadrics (SQ) [TUI IT3] . 

2. Simulating Hard Rigid Bodies with Sticky Spots 

HRBs can be simulated calculating their distance and collision contact 
point and time using Newton-Raphson method to solve a proper set of non 
linear equations as shown in [13]. In [13] we have also shown how to make 
use of linked lists (LL) and nearest neighbor lists (NNL) to attain good 
performance. If HRBs have a pronounced non-spherical shape, i.e. a large 
aspect ratio, it is mandatory to use NNL in order to minimize number of 
collision predictions. 

The event- calendar has been implemented using an hybrid approach as 
explained in [13] . where bounded priority queue is built on top of a binary 
tree implemented as in [15]. All operations (insertion, deletion and next 
event retrieve) with such priority queue have a complexity 0(1) with respect 
to particles number. In order to avoid round-off problems for events very 
close to each other, we shift forward time origin periodically. 

Particles interact through an hard core potential and they may be deco- 
rated with spherical sticky spots (SS) that interact via a square well potential 
[8] . Predicting next possible collision means that we have to take into account 
both collisions between spots and hard core collisions between HRBs. More 
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specifically if t$w is the next collision time between spots of two particles A 
and B and tnc is the next collision time between their hard cores, the next 
event between A and B will occur at time t next = mm{tsw, tuc}- All the 
details for an efficient algorithm to find collision times between spots can be 
found in [BJ. 

SSs have to be considered in calculating escape time from NNL. If t b ^ c 
is the escape time of the HRB from its bounding box (BB) and tg S is the 
escape time of its SSs, then the escape time tsB of a particle from its BB will 
be tBB = m in{^ifC' ^ssl- ^hc can De calculated evaluating the time evolution 
of the distance between the HRB and its BB, i.e. d(t) = minj{dj(t)}, where 
i = 1 ... 6 labels each of the 6 BB planes and di(t) is the distance between 
the given HRB and z-th plane. 

If the HRB is decorated with SSs we have also to consider the escape time 
of SSs from its surrounding BEQ The escape time can be calculated in this 
case considering simultaneous evolution of all distances between SSs and the 
6 BB planes. 

3. A Novel Method to Compute Escape Time 

In [13] calculation of the escape time tjj C of an HRB from its enclosing BB 
requires the evaluation of collision time between the HRB and its BB. This 
is computationally quite expensive although much faster than the prediction 
of the collision time between two HRBs. A possible trick is to consider a 
polyhedron enclosing the HRB and to evaluate the escape time from its BB 
(see Fig{T]). This polyhedron must be chosen to fit the HRB, i.e. given a 
certain shape (e.g. a parallelepiped) it should have the smallest possible 
size in order to enclose completely the HRB. Latter requirement will ensure 
that the escape time of such polyhedron tg C will be an underestimate of the 
escape time t b ^ c , i.e. t s ^ b c < t b ^ c . Evaluation of t s ^ b c consists in calculating the 
smallest escape time of all vertices of the polyhedron. At this point note that 
such vertices can be thought as spots of null diameter (NSS) meaning that the 
same algorithm used for calculating the escape time of SSs can be employed 
to evaluate the escape time tf^ b c . In Fig. [ljthe polyhedron is a parallelepiped 
that encloses the given cylindrical-like HRB. If shape of particles is more 
complicated than the one show in Fig. [T] simply a polyhedron with a shape 



1 Note that in this case BB must enclose both the HRB and the SSs 
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Figure 1: SQ enclosed in its bounding box (cyan parallelepiped) with SS of null diameter, 
represented here as finite size yellow spheres. 



that fits better the HRB shape can be used, i.e. there is no need to calculate 
analytically any jacobian as in the method proposed in [13] . 

4. Nearest Neighbour Lists with Null Spots: Performance Results 

In this Section we will test performance of the novel method described in 
Section [3] We consider here superquadrics (SQ), whose surface is defined as 
follows: 

/(*, y, z) = \x/a\ n + \y/b\ m + \z/c\ p -1 = (1) 

where the parameters n,m,p are real numbers and a, b, c are the SQ semi- 
axes. A monodisperse system of N = 512 SQs has been simulated with 
n — 8, m = p = 2 and with two equal semi-axes, i.e. b = c. Such SQs 
can be characterized by the elongation X = a/b and if elongation X < 1 
particles are called "oblate", while if X > 1 particles are called "prolate". 
For this test we have taken into account only prolate SQs, whose shape 
resembles that of a cylinder with smoothed edges (see Fig. [I]). In particular 
elongations X Q = 1.0,2.0,3.0 have been studied. The system of prolate 
SQs has been simulated in a cubic box of volume V with periodic boundary 
conditions at the volume fractions (ft = 0.20,0.30,0.35,0.40 (0 = irX b 3 p/6, 
where p = N/V is the number density). The length of the smallest semi-axes 
is chosen to be the unit of lenght (6 = 1.0), the mass of the SQ is the unit of 
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Figure 2: (a) Speedup Ssq versus elongation X for <j> = 0.20,0.40 (b) Speedup Ssq 
versus <j) for X = 1.0, 3.0. 



mass (m = 1) and the moment of inertia is spherically symmetric and equal 
to 1.0. To create the starting configuration at a desired 0, an extremely 
diluted crystal has been melted; afterwards, the particles have been grown 
independently up to the desired packing fraction (quench in at fixed N, 
Xq), similarly to what was done in [T3] . 

To test the algorithm speed, we use the CPU tim^j per collision, i.e. 
r c = T-tat/N con, where T tot is the (real) time needed to perform N co u collisions 
during a simulation. 

We can define the speedup Ssq as follows: 

Ssq = t c n /t? (2) 

where N refers to simulations that use NNL and S refers to simulations 
that make use of SNL. Fig. [2] (a) shows Ssq as a function of elongation 
Xq for two different volume fractions and it is apparent that use of SNL 
offers a speedup around 2 independently of elongation. In a similar way Fig. 
[2] (b) shows speedup Ssq for two elongations X = 1.0, 3.0 as a function 
of and it turns out again that Ssq is also independent of 0. Hence the 
novel method proposed for calculating the escape time of a HRB from its 
BB provides a significant speedup over previous methods and this speedup 
is quite independent of Xq and <fi. This result can be fairly understood 



2 CPU time means the real time spent by the CPU for calculations. 
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because number of collision predictions for a given X and (j> are independent 
of method used to calculate the escape time (either SNL or NNL). For that 
reason algorithm performance depends mainly on efficiency in rebuilding SNL 
or NNL, but SNL are more efficient than NNL and CPU time needed to 
calculate escape time is quite insensitive to elongation and volume fraction. 

5. Binary Mixtures of Hard-Sphers: Multiple Linked Lists 

Linked lists^] (LL) [15] are commonly employed in molecular dynamics 
simulations in order to avoid to check all the N 2 possible collisions among 
N simulated particles. In the LL method, the simulation box is partitioned 
into M 3 cells and only collisions between particles inside the same cell or 
its 26 adjacent cells are accounted for. This also means that, whenever an 
object crosses a cell boundary going from cell a to a new cell b, it has to 
be removed from cell a and added to cell b. Several variants of original LL 
method have been proposed in literature [16] with a view to improving the 
original LL method (as discussed for example in [17]). All these improved 
LL methods are intended to avoid unnecessary distance calculations but they 
do not tackle the case of a system composed of different species having very 
different sizes. Here we consider the case of a system composed of N s several 
species having different sizes {<Ti, . . . <Jn s } and for the sake of simplicity we 
can assume that particles are spheres. 

For making use of the original LL method, the simulation box must be 
partitioned into cubic cells and each cell must have a side length slightly 
greater than a max = maxj{<7j}. In general, this restriction can considerably 
compromise algorithm performance. To understand this, consider par- 
ticles of diameter Oi inside a cell and = (a ma x/&if , at a fixed volume 
fraction, if ^ — >- then — > oo and performance are severely compromised. 

A possible generalization of the original LL method for dealing with such 
polydisperse system consists in using one different LL for each pair of species, 
meaning that for each pair of species a different partitioning of simulation 
box into cells is used as illustrated pictorially in Fig. [3] for a binary mixture 
of hard spheres. 

In general if one has N s different species, N S (N S + l)/2 different LLs must 
be used where the cell side length for interacting species / and m F] has to be 



3 they are also referred to as linked cell lists 
4 Assuming that they are additive. 
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Figure 3: Pictorial representation of Multiple Linked Lists for a binary mixture of two 
particles where diameter a a of particles A is bigger than diameter ub of particles B. For 
each possible interaction AA, AB and BB a different partitioning of simulation box into 
cells is employed, i.e. multiple linked cell lists are used. 



greater than (07 + er m )/2. Implementation of MLL is quite straightforward 
except for the following remark: if particles, as they cross simulation box 
boundaries, are reinserted into the box with periodic boundary conditions, 
one has to ensure that reinsertion happens only once. 

6. Speedup of Multiple Linked Lists 

In this Section we will test performance of MLL method for a binary 
mixture of spheres having very different size. The two species (labeled by 
A and B) are characterized by a diameter ratio ua/^b = q > 1 and their 
masses are chosen to be equal and unitary. Na and Nb will be number of 
particles A and B respectively and iV = Na + Nb- The number of particles 
A will be kept fixed to 250, i.e. Na = 250. We will investigate the algorithm 
performance varying N, the volume fraction and the size ratio q. The 
simulation box is cubic with periodic boundary conditions. As for SQs we 
make use of the 0(1) event- calendar proposed in [14J. 

Again we define the speedup Ssq as follows: 

Sbm = r c LL /r c MLL (3) 

where LL refers to simulations that use single LL, MLL refers to simulations 
that employ MLL and r c is the CPU time per collision. Fig. [1] (a) shows 
Sbm as a function of 4> for two different size ratios q. It is remarkable that 
speedup ranges from 15 to 40. 
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Figure 4: (a) Speedup Sbm versus total volume fraction for binary mixtures for q — 2,5 
with N = 30000 for all points, (b) Speedup S B m versus size ratio q for N = 10000, 30000 
with <p — 0.60 for all points. Dashed lines are fits to function S l g M defined in Eq. |5|. 



The number of spheres of type B, Ng U , within a cell using conventional 
LL method is roughly: 



N cell = / 4 x 

vr [q%l - </>)N A /N B + 1] 1 1 
and a analytical estimate S^ M for the speedup turns out to be ^ 

S$ M = KN B el1 (5) 



where K is an arbitrary constant. Fig. [4] (a) and (b) shows also fits to 
simulations data of function S^ M defined in Eq. Agreement between 

numerical data and proves thctt Sj^m 

is a reasonable estimate. 



7. Conclusions 

In this paper two novel techniques easy to implement have been proposed 
for optimizing EDMD simulations. The SNL method offers a nearly con- 
stant speedup around 2 with respect to the old method proposed in [13] and 
it can be easily adapted to more complicated shapes of simulated particles. 
The SNL method is actually used for simulating a recently developed model 



5 for MLL method N% eU is of the order of 1. 
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of DNA duplexes (DNAD) [T5| consisting in cylindrical-like SQs decorated 
with two sticky spots on their two bases in order to model stacking interac- 
tions between DNADs. The details of latter study will be given in a future 
publication. 

MLLs have been already used in [Tj5] where spherical particles, evolv- 
ing according to event-driven brownian dynamics [201, can be absorbed by 
one fixed spherical sink whose size may be much greater than diffusing par- 
ticles. MLLs are also currently used for investigating phase diagram of a 
binary mixture of hard spheres whose size ratio q = 5 upon changing the 
partial volume fractions of two species. For this system it is possible to cal- 
culate a theoretical phase diagram with respect to glass transition within the 
framework of Mode Coupling Theory and recently jamming lines for such 
system have been evaluated both theoretically (using Replica theory) and 
numerically [2 lj. In view of growing interest for this system it is crucial to 
have an efficient algorithm to explore whole phase diagram. It is worth not- 
ing that also conventional time-driven molecular dynamics may benefit from 
MLL. Huge increase in performance provided by MLL, when the length scales 
between the components is so important, could play a relevant role in mul- 
tiscale simulations, a topic that is attracting a lot of interest for biological 
and material applications. 

Finally the two optimization methods, SNL and MLL, illustrated in this 
paper can be used together, like LL and NNL, by using MLL to generate 
SNL. 
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